According to the phase field theory, the governing equation for phase transition is written as 1 where τ₀ is the kinetic characteristic time, a_s(n) is the anisotropy function of the interfacial energy. n is the normal vector directing from solid to liquid. It is a function of space and time, i.e., n ≡ n(x, t). ϕ is a dimensionless order parameter, which is a continuous function depending on (x, t) and defined over the whole computational domain. W₀ is the thickness width of the diffusion-type interface. In the present work, we assume that ϕ = −1 and ϕ = + 1 correspond respectively to the liquid phase and the solid phase. V is an anisotropic vector related to a_s(n), and Q_ϕ represents the source term coupling with the heat transfer equation.

Through the Chapman–Enskog expansion, the phase field equation can be discretized into the LB scheme with anisotropic coefficients. The governing LB equation with discrete velocities for phase change is expressed as^[20,21] 2 where g_α(x, t) is the discrete distribution function for the phase field at space x and time t, and represents the equilibrium distribution function. τ_P (x,t) is the relaxation time for phase transition. It depends on the interfacial anisotropy function through . Because τ_P(x,t) is a function of space and time, we should update it in each time step on all lattices. It is obvious that the transport term of Eq. (1) is non-local and nonlinear at a time. This term is converted into local and linear operations as suggested in Eq. (2). In the the LB–PF models previously developed^[16–18] the transport term is solved by finite different method, while the governing equation can be solved by the linear (streaming) and local (collision) procedures in the present model. The linear and local features entitle the present model to performing parallel computation easily. Because the LB equation fits into the standard LB framework with the anisotropic coefficient , it is called the anisotropic LB–PF scheme. Therefore, the anisotropic LBE enables reasonable simulations of crystal growth with anisotropic patterns.

The equilibrium distribution function for the anisotropic LBE is written as 3 where v_n is the interface advancing velocity evoked by the surface energy. It is given as . Here, V is computed by 4 By computing the zero-order moment of distribution function, we obtain the order parameter ϕ at (x, t) through ϕ(x, t) = ∑_α g_α (x, t). Usually, the anisotropy function for a cubic crystal system is given as 5 Here, n is computed through n(x, t) = − ∇ ϕ(x, t)/|∇ ϕ(x, t)|. The present anisotropy function is used to simulate crystal growth with preferring orientation which is parallel to the x axis. Coordinate rotation should be made to simulate crystal growth with various preferring orientations. The last term of Eq. (2) is the source term which involves the thermal undercooling as the driven force to induce phase changes. The Q_P(x,t) in the source term is defined by 6 where λ represents the coupling coefficient, and is the dimensionless temperature. The is computed through with the specific heat c_p, the temperature T(x,t), the melting temperature T_m and the latent heat L_h.

The heat transport during phase change is driven by both melt flows and temperature gradients. Evolution of flow field is therefore modeled to provide velocities used in convective–diffusion equation for heat transport. Considering the multi-relaxation-time LB model has better numerical stability than the single-relaxation-time LB model, the general MRT–LB scheme is used to model melt flows with phase changes. The MRT–LB equation on a D-dimensional lattice is^[25,26] 7 where f_α is the distribution function representing the probability of finding a pseudo-particle of fluid in the α-th direction. is the corresponding equilibrium distribution function. e_α and δ t are discrete velocity and time step, respectively. M is a Q × Q orthogonal transformation matrix. Q represents all directions of the probable velocities on a D-dimensional lattice. S is a non-negative Q × Q diagonal matrix consisting of (s₀,s₁,…,s_{Q − 1}). Here, the subscripts α and β are used to represent components of vectors and matrices. A simplified form of transformation matrix M is used as suggested in Ref. [27]. The details for the calculation of can be found in our previous work^[23] and Ref. [25].

The D2Q9 model is used to discrete the two-dimensional space for its numerical stability and computational accuracy. The discrete velocities are expressed as 8 where the lattice speed c can be computed by c ≡ Δx/Δt. The lattice sound speed c_s is correlated to the lattice sound speed c by . The diagonal matrix S is expressed as 9 The relaxation rates, s_ν and s_e, are related to the kinematic shear viscosity ν and the bulk viscosity ζ, respectively. For the D2Q9 discrete velocities, the relationships between relaxation rates and viscosities are 10 The elements s_ε and s_q are flexible parameters which can be set freely in the range (0, 2). However, simulations may loss stability without optimized values. Therefore, linear analysis is usually used to obtain optimized parameters and achieve good numerical stabilities. In this work, the parameters are chosen as s_e = 1.19, s_ε = 1.4, s_q = 1.2 for all the simulations.

The multi-distribution LB approach is employed to simulate heat transfer during crystal growth. Another set of pseudo-particle distributions is used to describe the temperature field. The LB equation with the single relaxation time and a source term is written as 11 where h_α is the pseudo-particle distribution function, τ_T is the dimensionless relaxation time which is related with the thermal diffusion coefficient α through , and Q_T(x,t) represents the heat source caused by phase change. The source term is computed by Q_T(x, t) = ∂ϕ/2∂ t. The equilibrium distribution function depends on the local temperature and flow velocity u: 12 The dimensionless macroscopic temperature is calculated by summarizing the pseudo-particle distribution function: .

For a certain physical problem, it is usually given macroscopic variables at boundaries. However, the LB method uses microscopic particle distribution functions to model evolutions of macroscopic fields. Therefore, the particle distribution functions at boundary nodes should be computed through the proper numerical scheme according to the known conditions. Because the present model is a diffusion-type interface model, the bounce-back rule^[28,29] commonly used in sharp interface models is not available to impose the non-slip boundary condition. In the framework of LB method, there are several boundary conditions to treat the S–L interface. When considering diffusion interfaces, the porous medium no-slip boundary condition is usually incorporated into the LB scheme. We therefore treats the diffusive solid–liquid interface as a kind of porous medium by the following scheme 13 where represents the increment of α-th distribution function due to relaxation, and accounts for the increment of α-th distribution function by bounce-back. The is computed by 14 where φ_s is the local solid phase fraction, and represents the opposite direction of the α-th direction. The solid fraction is calculated by ϕ by φ_s = (ϕ + 1)/2.

It has been proofed that the present model is equivalent to the model developed by Beckermann^[7] for simulations of crystal growth with melt convection. Details of model validation involving simulations of free dendritic growth into an undercooling melt in pure diffusive and convective conditions can be found in our previous work.^[22]

Pattern evolution during crystal growth from undercooling melts was firstly investigated in the condition of pure conductive condition. In the simulation, a circular crystal seed was initially placed at the center of an L × L domain. The growth preferring orientation of the seed is aligned with the x axis. The periodic boundary conditions were imposed on the four sides of the temperature field, i.e., and . The mesh of 2048 × 2048 was used in the simulation. The spacial step and time intervals were set as δ x = 0.4W₀ and δ t = 0.008 τ₀, respectively. The capillary length d₀ = 0.185 W₀ and the thermal diffusion coefficient with the coupling coefficient λ = a₁ W₀/d₀, where a₁ = 0.8839.

Figure 1 shows the simulated crystal patterns at time t = 512 τ₀ with d₀ = 0.185W₀, 15ε = 0.15, and various initial undercoolings. Contours present the the solid–liquid interface with ϕ = 0. As shown in the figure, the temperatures at center of the domain are higher than the boundaries because the latent heat was released during the crystal growth. Despite different initial undercoolings, the temperature fields and the interface contours display four-fold symmetrical patterns. In the case of , the shape of the growing crystal pattern was a dendritic morphology. When the initial undercooling was increased, the crystal kept growing as a dendrite but the doublons generated during the growing process, as shown in the case of . With the increase of the initial undercooling, the seaweed morphology can be reproduced numerically. For the case , the increased undercooling significantly modified the growth pattern. Each single primary arm was split into triple arms. When the initial undercooling was increased to a higher value, , the growing crystal demonstrates cauliflower-like pattern.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. Simulated crystal morphologies in the pure conductive condition at time (a) t = 2560τ0, (b) t = 1280τ0, (c) t = 640τ0, (d) t = 320τ0 with d0 = 0.185, ε = 0.01 and various initial temperatures (a) , (b) , (c) , (d) .

In order to quantitatively characterize the crystal pattern, the angle θ_n, which is defined as the angle between the directions of the interface normal vector and the x axis, is measured for the cases of various initial undercoolings. The range of θ_n is [−180°, 180°]. Figure 2 depicts the distributions of interface angle θ_n at time t = 512τ₀ with d₀ = 0.185W₀, ε = 0.01, and various initial undercoolings. As shown in Fig. 2, the angles demonstrate periodic distributions, which reflects the symmetrical pattern of crystals. For the case of , the crystal patterns obviously present four-fold symmetrical features, and there are four main-valleys and four main-peaks of in the frequency-angle histogram. When , four doublons appear, and the secondary valleys and peaks correspondingly become more clearly in the histogram. With the increase of initial undercooling, the growing crystal transforms from dendritic pattern into seaweed and cauliflower-like patterns, which result in a more complicated S–L interface. Correspondingly, the periodicity of angle distribution becomes less obviously, as shown in Figs. 2(c) and 2(d). In spite of this, the tiny valleys and peaks can also be observed from the F_θ–θ_n histograms. For the case of , the small secondary deepest valleys indicate the secondary arms, as shown in Fig. 1(c) and Fig. 2(c). For the case of , no obvious arm can be found from the crystal pattern. Correspondingly, the peaks and valleys become less distinctly, as shown in Fig. 2(d). The standard deviations of the four θ_n-distributions in Fig. 2 are 0.0235, 0.0259, 0.0248, and 0.0135. A larger standard deviation indicates that the crystal grew with obvious preferred orientations, and the crystal would develop into a dendritic morphology. On the contrary, a smaller standard deviation suggests the seaweed morphology. Therefore, the distribution of interface normal angles reflects the symmetrical feature of crystal morphologies. In all the simulated cases, a larger standard deviation of θ_n-distribution corresponds to the dendritic morphology, while a smaller deviation of θ_n-distribution provide the evidence of seaweed patterns.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. Interface normal angle distribution for crystal growth in the pure conductive condition at time (a) t = 2560τ0, (b) t = 1280τ0, (c) t = 640τ0, (d) t = 320τ0 with d0 = 0.185, ε = 0.01 and various initial temperatures (a) , (b) , (c) , (d) .

Figure 3 shows the interface curvature distribution for the crystal growth cases in Fig. 1. The K-distribution demonstrates high and narrow shapes at a relative small undercooling, which indicates the dendritic and seaweed morphologies. While at the largest undercooling, the K-distribution exhibits broad and smooth reflecting the cauliflower-like pattern. The standard deviations of the K-distribution in Figs. 3(a)–3(d) are 0.02353, 0.02592, 0.02480, and 0.01355, respectively. The first three deviations are very close to each other, and there values are much higher than the last one. Therefore, increasing the undercooling would change the interface curvature distributions and result in the decrease trends of the K-distribution deviations. With the increase of the undercooling, there is a evolution of crystal morphology from dendrite, seaweed to cauliflower. In all the cases, the narrowest distribution corresponds to the dendritic morphology while the broadest distribution reflects the cauliflower-like pattern.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. Interface curvature distribution for crystal growth in the pure conductive condition at time (a) t = 2560τ0, (b) t = 1280τ0, (c) t = 640τ0, (d) t = 320τ0 with d0 = 0.185, ε = 0.01, and various initial temperatures (a) , (b) , (c) , (d) .

To analyze the convective effect of melt flows on pattern evolution, we simulated the crystal growth in conditions of d₀ = 0.185, ε = 0.01, and various initial temperatures. The left side of the domain is supposed to be inlet and the right side is outlet. At the inlet the velocity is given as u_in = 1.00, and at the outlet the boundary condition is ∂ u/∂ x = 0. Here, u_in = 1.00 is the incoming flow velocity at the inlet of the computational domain, whose unit is W₀/τ₀. We ignored the velocity unit to simplify expressions in the following context. The non-equilibrium extraordinary scheme is used to deal with the inlet and outlet boundary conditions of the flow field. Figure 4 displays the simulated crystal morphologies under several initial undercoolings, , −0.70, −0.80, and −0.90. Other conditions were set the same as those in Fig. 1.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. Simulated crystal morphologies at the total solid fraction in the condition of melt flows with d0 = 0.185, ε = 0.01, uin = 1.00, and various initial temperatures (a) , (b) , (c) , (d) .

Compared the morphologies in Fig. 4 with Fig. 1, it can be found that the melt flows can significantly change the crystal growth and final patterns. As shown in Fig. 4, the melt flows promote the growth of upstream arms and retarded the growth of downstream arms. The upstream crystal arms grew more quickly and became longer and thicker than the arms in downstream regions. With the impingement of melt flows, the temperature gradient around the upstream tip is larger than those around the vertical and downstream tips. The temperature fields are thus significantly distorted by the melt flows. As a result, the upstream and downstream arms become asymmetrical, which is different from the no-flow condition. Because the flow field is symmetrical, and the temperature field keeps a symmetrical distribution. Therefore, the driven force at the S–L interface is symmetrical about the line y = 1024δ x, and the final crystal morphologies keep axial symmetrical shapes despite of different initial undercoolings.

It can also be found that the vertical arms slightly display the upstream growing trend. The growing orientations of vertical arms are not exactly the vertical direction but the direction tilted to upstream. That is because the heat fluxes around vertical arms on the upstream side are higher than on the downstream side. The growth speed at the upstream side is therefore higher than the downstream side. As a result, the vertical arms are growing against the incoming flows, and the vertical arms appear tilted by the melt flows. Interestingly, the tip-splitting can be captured by the present simulation, as shown in Figs. 4(c) and 4(c). When increasing the initial undercoolings, the crystal grows at a relative fast solidification rate. More latent heat is released at the interface in the condition of larger initial undercoolings. When the solidification rate reach to a critical value, the tip-splitting occurs and the crystal tip grows without following the crystallographic orientations. There exist an interface instability during crystal growth. Such instability can break up the growth symmetrical feature under certain conditions. When the crystal reaches a critical growing speed, The tip-splitting occurs when the interface loses stability.

Figure 5 depicts the interface normal angle distributions for the cases in Fig. 4. All the histograms are symmetrical about the line θ_n = 0°, which corresponds to the vertical symmetrical features of the crystal morphologies in Fig. 4. At the higher initial temperatures, such as and −0.70, four deep valleys can be clearly observed around θ_n = 0°, ± 90° and ± 180°, which means that the crystal arms are well developed along the directions. When the initial temperatures is increased, the histogram peaks and valleys become lower. It indicates the transition of crystal morphologies from dendrites, seaweed to cauliflower-like patterns, which is the same as that in the no-flow conditions. Interestingly, the upstream growing trend is quantitatively captured by the histograms. For example, the slope |dF_θ/dθ_n| at the angles θ_n = 90° + δ θ_n is larger than the slope at the angles θ_n = 90° − δ θ_n. Here, δ θ_n represent a tiny angle value. Figure 6 shows the histograms of interface curvature distributions. Compared the K-distributions with the cases in the pure conductive conditions, we can find that the melt flow does not change the pattern transition when decreasing initial temperatures. The histograms become lower and broader with the increase of initial undercoolings, which reflects the morphological transition as shown in Fig. 1.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. Interface normal angle distribution for crystal growth at the total solid fraction in the condition of d0 = 0.185, ε = 0.01, uin = 1.00, and various initial temperatures (a) , (b) , (c) , (d) .

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. Interface curvature distribution for crystal growth at the total solid fraction in the condition of Fs = 15.5% with d0 = 0.185, ε = 0.01, uin = 1.00, and various initial temperatures (a) , (b) , (c) , (d) .

Figure 7 shows the simulated morphologies with ε = 0.01, , and various incoming flow velocities. The incoming flow field was imposed the inlet velocities of u_in = 0.05, 0.50, 1.00, and 2.00 in units of W₀/τ₀. With the increase of flow strength, more and more undercooled melt and released latent heat were flushed away from the upstream to downstream regions, and the upstream arms have larger driven forces to grow. Subsequently, the asymmetrical feature in the horizontal direction became more and more obvious. We can also find that the vertical arms display the trends of growing upstream-wards when u_in = 1.00 and 2.00, as shown in Figs. 7(c) and 7(d). As the upstream-wards growing trend is caused by the melt flow which tilts the direction of heat fluxes to upstream direction, the trend therefore becomes more and more obvious with the increase of the flow strength. The tip-splitting can also be found in the simulations of higher inlet velocities. As discussed above, the tip-splitting is evoked by the interface instability when the solidification rate reach a critical value. The higher of the solidification rate, the more instability the interface will have. The crystal growing from melts of higher inlet velocities trends to form morphologies with more extensively split tips.

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. Simulated morphologies for crystal growth at time t = 512τ0 with ε = 0.01, , and various incoming flow velocities (a) uin = 0.25, (b) uin = 0.50, (c) uin = 1.00, (d) uin = 2.00.

Figure 8 shows the interface normal angle distributions corresponding to the cases in Fig. 7. By comparing with the simulated morphologies in Fig. 7, it can be found that the peaks and valleys of the θ_n-distributions in Fig. 8 reflects the crystal growth orientations and doublons. For example, the valleys V₁ and V₂ in Fig. 8(c) indicate the arms A₁ and A₂ as show in Fig. 7(c). With the incoming flow velocity being increased, the peaks and valleys appear more clearly reflecting to the crystal more obviously transverse-asymmetrical morphologies. Despite of various melt flow strength, all the θ_n-distributions keep symmetrical respecting to the line of θ_n = 0. It suggests the crystal morphologies are symmetrical in the longitudinal direction.

Fig. 8.

	Figure Option ViewDownloadNew Window
	Fig. 8. Interface normal angle distribution for crystal growth at time t = 512τ0 with ε = 0.01, , and various incoming flow velocities (a) uin = 0.25, (b) uin = 0.50, (c) uin = 1.00, (d) uin = 2.00.

Figure 9 shows the interface curvature distributions for the crystal growth in Fig. 7. In all the cases, the K-distribution histograms display high and narrow distributions, from which one can deduce the crystal would grow up into a seaweed morphology. In spite of different incoming flow velocity, the curvature distributes in the region around K = 0. It means that there are a lot of flat interfacial area during the seaweed growth. Figure 8 provides the evidence of this deduction.

Fig. 9.

	Figure Option ViewDownloadNew Window
	Fig. 9. Interface curvature distribution for crystal growth at t = 512τ0 with ε = 0.01, , and various incoming flow velocities (a) uin = 0.25, (b) uin = 0.50, (c) uin = 1.00, (d) uin = 2.00.

Table 1.

Table 1.

Table 1. The mean value FK,m and standard deviation SK of the curvature distributions for crystal growth at t = 512τ0 with ε = 0.01, , and various flow velocities uin. . Flow velocity, uin Mean value, FK,m Standard deviation, SK 0.00 0.01000 0.02235 0.25 0.01000 0.02064 0.50 0.00999 0.02091 1.00 0.00983 0.01682 2.00 0.00951 0.01712

Table 1. The mean value FK,m and standard deviation SK of the curvature distributions for crystal growth at t = 512τ0 with ε = 0.01, , and various flow velocities uin. .

The peaks of K-distributions were changed with the increase of the flow strength. Table 1 lists the mean value and standard deviation of the K-distributions for all the cases with different flow velocity u_in. It can be found that the mean value and the standard deviation of curvatures slightly decrease when enhancing the flow strength. The curvature distribution has a higher dispersivity in the condition of a larger u_in. That is because that the melt flow changed the temperature field significantly, and the crystal would grow into more complex morphologies with the increase of flow strength.

In addition, we carried out simulations of crystal growth with ε = 0.0 to examine the usage of θ_n- and K-distributions characterizing morphologies. Figure 10 displays the simulated morphologies. It can be seen that at the lowest flow strength the crystal grows up with a smooth interface. Neither split nor branch can be observed clearly from the morphology. When increasing the incoming flow velocity to u_in = 0.25, the flow field evokes the instability of the advancing interface. When u_in = 0.50, the two splits can be clearly observed, as shown in Fig. 10(c). In the condition of the largest flow velocity, the flow field raises the strongest interface instability. As a result, splits and branches generate at the interface in the upstream region. Correspondingly, as shown in Fig. 11, peaks and valleys appear in the interface normal angle distribution histograms. Even though the flow strength is different for all the cases, the symmetrical feature of the crystal growth can be observed both from the figures of crystal morphologies and θ_n-distribution histograms.

Fig. 10.

	Figure Option ViewDownloadNew Window
	Fig. 10. Simulated morphologies for crystal growth at time t = 512τ0 with ε = 0.0, , and various incoming flow velocities , (a) uin = 0.25, (b) uin = 0.50, (c) uin = 1.00, and (d) uin = 2.00.

Fig. 11.

	Figure Option ViewDownloadNew Window
	Fig. 11. Interface normal angle distribution for crystal growth at time t = 512τ0 with ε = 0.0, , and various incoming flow velocities (a) uin = 0.25, (b) uin = 0.50, (c) uin = 1.00, (d) uin = 2.00.

Figures 12 displays the interface curvature distribution for the crystal growth. As shown in the figure, the curvature distributes narrowly and highly at the lower flow strength. At the higher flow strength, the crystal started to generate branches and splits, and the curvature distributes dispersively. The K-distribution trends of enlarging flow velocity in the present condition is the same as that in the condition of ε = 0.01. Therefore, the θ- and K-distributions reflect crystal morphological features and can be used to describe quantitatively the evolution of crystal patterns.

Fig. 12.

	Figure Option ViewDownloadNew Window
	Fig. 12. Interface curvature distribution for the crystal growth at time t = 512τ0 with ε = 0.0, , and various incoming flow velocities (a) uin = 0.25, (b) uin = 0.50, (c) uin = 1.00, (d) uin = 2.00.